!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000, 2003 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!###########################################################################

SUBROUTINE radiate(mzp,mxp,myp,ia,iz,ja,jz,mynum)

  USE mem_tend   ,  ONLY: tend
  USE mem_grid   ,  ONLY: ngrid,time,dtlt,itime1,nzg,nzs,npatch,grid_g,nnzp, &
                          if_adap,zm,zt
  USE mem_leaf   ,  ONLY: leaf_g
  USE mem_radiate,  ONLY: ilwrtyp,iswrtyp,radiate_g,radfrq, &
                           ncall_i,prsnz,prsnzp
  USE mem_basic  ,  ONLY: basic_g
  USE mem_scratch,  ONLY: vctr1,vctr2,vctr3,vctr4,vctr5,vctr6,vctr7,vctr8, &
                          vctr9,vctr10,vctr11,vctr12,scratch
  USE mem_micro  ,  ONLY: micro_g
  USE rconstants ,  ONLY: cp,cpor,p00
  USE rrad3      ,  ONLY: jday,solfac,ng,nb,nsolb,npsb,nuum,prf,alpha,trf,beta,&
                          xp,wght,wlenlo,wlenhi,solar0,ralcs,a0,a1,a2,a3, &
		 	  exptabc,ulim,npartob,npartg,ncog,ncb,ocoef,bcoef,gcoef
  USE ref_sounding, ONLY: pi01dn

  ! CATT
  !kmlnew
  USE micphys     , ONLY: gnu,level,icloud,irain,ipris,isnow,iaggr,igraup,ihail
  USE mem_cuparm  , ONLY: cuparm_g, nnqparm
  !kmlnew
  USE rad_carma   , ONLY: radcomp_carma
  USE catt_start  , ONLY: catt           ! intent(in)
  USE mem_scalar  , ONLY: scalar_g

  IMPLICIT NONE

  INTEGER :: mzp,mxp,myp,ia,iz,ja,jz,koff,mynum

  REAL, pointer :: rc_ptr

  REAL :: solc

  !kmlnew
  REAL, DIMENSION(mzp,mxp,myp) :: LWL,IWL,fracao_liq
  REAL, DIMENSION(mxp,myp)     :: RAIN
  REAL :: dummy
  INTEGER :: ncall = 0
  INTEGER :: i,j,k
  !kmlnew
  INTEGER :: t1,c_rate,t2,ti
  DOUBLE PRECISION :: tg


!  CALL system_clock(ti, c_rate)


  IF (ilwrtyp + iswrtyp .eq. 0) RETURN

  CALL tEND_accum(mzp,mxp,myp,ia,iz,ja,jz,tEND%tht(1)  &
       ,radiate_g(ngrid)%fthrd(1,1,1))

!!$  PRINT *, "DEBUG: radiate: tEND_accum FEITO!"
!!$  CALL FLUSH(6)

  IF (mod(time + .001,radfrq) .lt. dtlt .or. time .lt. 0.001) THEN

!     PRINT 90,time,time/3600.+(itime1/100+mod(itime1,100)/60.)
!90   format('   radiation tendencies updated    time =',f10.1,  &
!          '  UTC TIME (HRS) =',F6.1)
	
!!$     CALL FLUSH(6)

     ! Compute solar zenith angle, multiplier for solar constant, sfc albeDO,
     ! and surface upward longwave radiation.

     !PRINT*, 'at a1 ',ngrid

     CALL radprep(mxp,myp,nzg,nzs,npatch,ia,iz,ja,jz,jday   &

     ,leaf_g(ngrid)%soil_water      (1,1,1,1)  &
          ,leaf_g(ngrid)%soil_energy     (1,1,1,1)  &
          ,leaf_g(ngrid)%soil_text       (1,1,1,1)  &
          ,leaf_g(ngrid)%sfcwater_energy (1,1,1,1)  &
          ,leaf_g(ngrid)%sfcwater_depth  (1,1,1,1)  &
          ,leaf_g(ngrid)%leaf_class      (1,1,1)    &
          ,leaf_g(ngrid)%veg_fracarea    (1,1,1)    &
          ,leaf_g(ngrid)%veg_height      (1,1,1)    &
          ,leaf_g(ngrid)%veg_albeDO      (1,1,1)    &
          ,leaf_g(ngrid)%patch_area      (1,1,1)    &
          ,leaf_g(ngrid)%sfcwater_nlev   (1,1,1)    &
          ,leaf_g(ngrid)%veg_temp        (1,1,1)    &
          ,leaf_g(ngrid)%can_temp        (1,1,1)    &

          ,solfac  &
          ,grid_g(ngrid)%glat       (1,1)  &
          ,grid_g(ngrid)%glon       (1,1)  &
          ,radiate_g(ngrid)%rshort  (1,1)  &
          ,radiate_g(ngrid)%rlong   (1,1)  &
          ,radiate_g(ngrid)%rlongup (1,1)  &
          ,radiate_g(ngrid)%albedt  (1,1)  &
          ,radiate_g(ngrid)%cosz    (1,1)  )


     CALL azero(mzp*mxp*myp,radiate_g(ngrid)%fthrd(1,1,1))

        ! Carma radiation
     IF (ilwrtyp==4 .or. iswrtyp==4) THEN

           !kmlnew
          ! IF (nCALL == 0) THEN
          !    nCALL = 1
	  !srf- problemas com RAIN quando conveccao esta OFF
              CALL azero(mzp*mxp*myp,LWL)
              CALL azero(mzp*mxp*myp,IWL)
              CALL azero(mzp*mxp*myp,fracao_liq)
              CALL azero(mxp*myp,RAIN)
          ! ENDIF

           IF (level == 2) THEN
              call atob (mxp * myp * mzp,micro_g(ngrid)%rcp,LWL)
	
              IF (nnqparm(ngrid) > 0) THEN
                 DO j=ja,jz
                    DO i=ia,iz
                       ! conv  precip at mm/h
                       RAIN(i,j)= cuparm_g(ngrid)%conprr(i,j)* 3600.
                    ENDDO
                 ENDDO
              ENDIF
	
           ELSEIF (level >= 3) THEN

              IF (nnqparm(ngrid) > 0) THEN
                 DO j=ja,jz
                    DO i=ia,iz
                        ! conv + explic  precip at mm/h
                       RAIN(i,j)= (cuparm_g(ngrid)%conprr(i,j) + &
                            micro_g(ngrid)%pcpg(i,j)) * 3600.
                    ENDDO
                 ENDDO
              ENDIF

              CALL azero2(mzp*mxp*myp,LWL,IWL)
              IF(icloud>0) &
                   call accum (mxp * myp * mzp,LWL,micro_g(ngrid)%rcp)
	
              !kml	   IF(irain>0) 	call accum (mxp * myp * mzp,LWL,micro_g(ngrid)%rrp(1,1,1))

              IF(igraup>0) THEN
	
                 DO j=ja,jz
                    DO i=ia,iz
                       DO k=1,mzp
                          call qtc(micro_g(ngrid)%Q6(k,i,j), &
                               dummy,fracao_liq(k,i,j))
                       ENDDO
                    ENDDO
                 ENDDO
                 call ae1t1p1(mxp * myp * mzp,LWL,fracao_liq,  &
                      micro_g(ngrid)%rgp,LWL)
	
                 fracao_liq(:,:,:) = 1. - fracao_liq(:,:,:)

                 call ae1t1p1(mxp * myp * mzp,IWL,fracao_liq,&
                      micro_g(ngrid)%rgp,IWL)
	
              ENDIF

              IF(ihail>0) THEN
	
                 DO j=ja,jz
                    DO i=ia,iz
                      DO k=1,mzp
                          call qtc(micro_g(ngrid)%Q7(k,i,j),dummy,  &
                               fracao_liq(k,i,j))
                       ENDDO
                    ENDDO
                 ENDDO
                 call ae1t1p1(mxp*myp*mzp,LWL,fracao_liq, &
                      micro_g(ngrid)%rhp,LWL)
	
                 fracao_liq(:,:,:) = 1. - fracao_liq(:,:,:)

                 call ae1t1p1(mxp*myp*mzp,IWL,  &
                      fracao_liq,micro_g(ngrid)%rhp,IWL)
	
              ENDIF
	
              IF(iaggr>0) call accum (mxp*myp*mzp,IWL,micro_g(ngrid)%rap)

              IF(isnow>0) call accum (mxp*myp*mzp,IWL,micro_g(ngrid)%rsp)

              IF(ipris>0) call accum (mxp*myp*mzp,IWL,micro_g(ngrid)%rpp)
	
           ENDIF

           CALL radcomp_carma(mzp,mxp,myp,ia,iz,ja,jz,solfac  &
                ,basic_g(ngrid)%theta       &
                ,basic_g(ngrid)%pi0         &
                ,basic_g(ngrid)%pp          &
                ,basic_g(ngrid)%rv          &
		,RAIN,LWL,IWL               &
                ,basic_g(ngrid)%dn0         &
                ,basic_g(ngrid)%rtp         &
                ,radiate_g(ngrid)%fthrd     &
                ,grid_g(ngrid)%rtgt         &
                ,grid_g(ngrid)%f13t         &
                ,grid_g(ngrid)%f23t         &
                ,grid_g(ngrid)%glat         &
                ,grid_g(ngrid)%glon         &
                ,radiate_g(ngrid)%rshort    &
                ,radiate_g(ngrid)%rlong     &
                ,radiate_g(ngrid)%albedt    &
                ,radiate_g(ngrid)%cosz      &
                ,radiate_g(ngrid)%rlongup   &
                ,mynum                      &
                ,grid_g(ngrid)%fmapt        &
		,leaf_g(ngrid)%patch_area   &
		,npatch                     &
		)

     END IF

     IF (ilwrtyp .le. 2 .or. iswrtyp .le. 2) THEN

        ! IF using Mahrer-Pielke and/or Chen-Cotton radiation, CALL radcomp.

        CALL radcomp(mzp,mxp,myp,ia,iz,ja,jz,solfac  &
             ,basic_g(ngrid)%theta     (1,1,1)  &
             ,basic_g(ngrid)%pi0       (1,1,1)  &
             ,basic_g(ngrid)%pp        (1,1,1)  &
             ,basic_g(ngrid)%rv        (1,1,1)  &
             ,basic_g(ngrid)%dn0       (1,1,1)  &
             ,basic_g(ngrid)%rtp       (1,1,1)  &
             ,radiate_g(ngrid)%fthrd   (1,1,1)  &
             ,grid_g(ngrid)%rtgt       (1,1)    &
             ,grid_g(ngrid)%f13t       (1,1)    &
             ,grid_g(ngrid)%f23t       (1,1)    &
             ,grid_g(ngrid)%glat       (1,1)    &
             ,grid_g(ngrid)%glon       (1,1)    &
             ,radiate_g(ngrid)%rshort  (1,1)    &
             ,radiate_g(ngrid)%rlong   (1,1)    &
             ,radiate_g(ngrid)%albedt  (1,1)    &
             ,radiate_g(ngrid)%cosz    (1,1)    &
             ,radiate_g(ngrid)%rlongup (1,1)    &
             ,mynum)

     END IF

     IF (iswrtyp .eq. 3 .or. ilwrtyp .eq. 3) THEN

        ! Using Harrington radiation


        !PRINT*, 'at a4 ',ngrid

        IF (nCALL_i .eq. 0) THEN

           ! IF first CALL for this node, initialize several quantities & Mclatchy
           ! sounding DATA.


           !PRINT*, 'at a5 ',ngrid

           CALL radinit(ng,nb,nsolb,npsb,nuum,prf,alpha,trf,beta  &
                ,xp,wght,wlenlo,wlenhi,solar0,ralcs,a0,a1,a2,a3,solc  &
                ,exptabc,ulim,npartob,npartg,ncog,ncb  &
                ,ocoef,bcoef,gcoef,gnu)

           prsnz  = (pi01dn(nnzp(1)-1,1) / cp) ** cpor * p00
           prsnzp = (pi01dn(nnzp(1)  ,1) / cp) ** cpor * p00


           !PRINT*, 'at a6 ',ngrid

           CALL mclatchy(1,mzp,IF_adap,koff  &
                ,prsnz,prsnzp  &
                ,grid_g(ngrid)%glat       (1,1)  &
                ,grid_g(ngrid)%rtgt       (1,1)  &
                ,grid_g(ngrid)%topt       (1,1)  &
                ,radiate_g(ngrid)%rlongup (1,1)  &
                ,zm,zt,vctr1,vctr2,vctr3,vctr4,vctr5,vctr6,vctr7  &
                ,vctr8,vctr9,vctr10,vctr11,vctr12)

           nCALL_i = nCALL_i + 1
        END IF

        ! For any CALL, interpolate the mclatchy sounding DATA by latitude and
        ! season.


        !PRINT*, 'at a7 ',ngrid

        CALL mclatchy(2,mzp,IF_adap,koff  &
             ,prsnz,prsnzp  &
             ,grid_g(ngrid)%glat       (1,1)  &
             ,grid_g(ngrid)%rtgt       (1,1)  &
             ,grid_g(ngrid)%topt       (1,1)  &
             ,radiate_g(ngrid)%rlongup (1,1)  &
             ,zm,zt,vctr1,vctr2,vctr3,vctr4,vctr5,vctr6,vctr7  &
             ,vctr8,vctr9,vctr10,vctr11,vctr12)

        ! IF using Harrington radiation with moisture complexity LEVEL < 3,
        ! CALL radcomp3 which is a substitute driving structure to the bulk
        ! microphysics.

        IF (level <= 2) THEN
           IF (level == 2) THEN
              rc_ptr => micro_g(ngrid)%rcp(1,1,1)
           ELSE
              CALL azero(mzp*mxp*myp,scratch%vt3dp(1))
              rc_ptr => scratch%vt3dp(1)
           END IF


           !PRINT*, 'at a8 ',ngrid

           CALL radcomp3(mzp,mxp,myp,ia,iz,ja,jz  &
                ,grid_g(ngrid)%lpw        (1,1)    &
                ,grid_g(ngrid)%glat       (1,1)    &
                ,grid_g(ngrid)%rtgt       (1,1)    &
                ,grid_g(ngrid)%topt       (1,1)    &
                ,radiate_g(ngrid)%albedt  (1,1)    &
                ,radiate_g(ngrid)%cosz    (1,1)    &
                ,radiate_g(ngrid)%rlongup (1,1)    &
                ,radiate_g(ngrid)%rshort  (1,1)    &
                ,radiate_g(ngrid)%rlong   (1,1)    &
                ,basic_g(ngrid)%rv        (1,1,1)  &
                ,basic_g(ngrid)%dn0       (1,1,1)  &
                ,radiate_g(ngrid)%fthrd   (1,1,1)  &
                ,basic_g(ngrid)%pi0       (1,1,1)  &
                ,basic_g(ngrid)%pp        (1,1,1)  &
                ,basic_g(ngrid)%theta     (1,1,1)  &
                ,rc_ptr  )
        END IF

     END IF
  END IF

  !PRINT*, 'at a9 ',ngrid


  RETURN
END SUBROUTINE radiate

!*****************************************************************************

SUBROUTINE tend_accum(m1,m2,m3,ia,iz,ja,jz,at,at2)

  IMPLICIT NONE
  INTEGER :: m1,m2,m3,ia,iz,ja,jz,i,j,k
  REAL, DIMENSION(m1,m2,m3) :: at,at2

!!$  PRINT *, "DEBUG: tEND_accum"
!!$  PRINT *, "                  ja,jz,ia,iz,m1=", ja,jz,ia,iz,m1
!!$  CALL FLUSH(6)

  DO j = ja,jz
     DO i = ia,iz
        DO k = 1,m1

!!$           PRINT *,"j,i,k=", j, i, k
!!$           CALL FLUSH(6)

           at(k,i,j) = at(k,i,j) + at2(k,i,j)
        END DO
     END DO
  END DO

!!$  PRINT *, "DEBUG: tEND_accum: FIM!"
!!$  CALL FLUSH(6)

  RETURN
END SUBROUTINE tend_accum

!*****************************************************************************

SUBROUTINE radprep(m2,m3,mzg,mzs,np,ia,iz,ja,jz,jday,soil_water,soil_energy, &
                   soil_text,sfcwater_energy,sfcwater_depth,leaf_class ,     &
                   veg_fracarea,veg_height,veg_albedo,patch_area, &
		   sfcwater_nlev,veg_temp,can_temp,solfac,glat,glon,rshort, &
		   rlong,rlongup,albedt,cosz)

  IMPLICIT NONE

  INTEGER :: m2,m3,mzg,mzs,np,ia,iz,ja,jz,jday
  REAL :: solfac
  REAL, DIMENSION(m2,m3) :: glat,glon,rshort,rlong,rlongup,albedt,cosz
  REAL, DIMENSION(mzg,m2,m3,np) :: soil_water,soil_energy,soil_text
  REAL, DIMENSION(mzs,m2,m3,np) :: sfcwater_energy,sfcwater_depth
  REAL, DIMENSION(m2,m3,np) :: leaf_class,veg_fracarea,veg_height,veg_albeDO  &
       ,patch_area,sfcwater_nlev,veg_temp,can_temp

  INTEGER :: ip,i,j
  ! REAL :: c1,c2

  ! Compute solar zenith angle [cosz(i,j)] & solar constant factr [solfac].

  CALL zen(m2,m3,ia,iz,ja,jz,jday,glat,glon,cosz,solfac)

  ! Compute patch-averaged surface albeDO [albedt(i,j)] and up longwave
  ! radiative flux [rlongup(i,j)].

  CALL azero2(m2*m3,albedt,rlongup)
  DO ip = 1,np
     DO j = 1,jz
        DO i = 1,iz

           !PRINT*, 'at b1 ',i,j,ip,ngrid

           CALL sfcrad(mzg,mzs,ip               &
                ,soil_energy    (1,i,j,ip) ,soil_water      (1,i,j,ip)  &
                ,soil_text      (1,i,j,ip) ,sfcwater_energy (1,i,j,ip)  &
                ,sfcwater_depth (1,i,j,ip) ,patch_area      (i,j,ip)    &
                ,can_temp       (i,j,ip)   ,veg_temp        (i,j,ip)    &
                ,leaf_class     (i,j,ip)   ,veg_height      (i,j,ip)    &
                ,veg_fracarea   (i,j,ip)   ,veg_albeDO      (i,j,ip)    &
                ,sfcwater_nlev  (i,j,ip)                                &
                ,rshort         (i,j)      ,rlong           (i,j)       &
                ,albedt         (i,j)      ,rlongup         (i,j)       &
                ,cosz           (i,j)                                   )

        END DO
     END DO
  END DO
  RETURN
END SUBROUTINE radprep

!*****************************************************************************

SUBROUTINE radcomp(m1,m2,m3,ia,iz,ja,jz,solfac  &
     ,theta,pi0,pp,rv,dn0,rtp,fthrd  &
     ,rtgt,f13t,f23t,glat,glon,rshort,rlong,albedt,cosz,rlongup  &
     ,mynum)

!srf-rasm60  USE mem_grid   , ONLY : dzm,dzt,nzp,itopo,plonn,ngrid,time,itime1,centlon
  USE mem_grid   , ONLY : dzm,dzt,nzp,itopo,plonn,ngrid,time,itime1,centlon,grid_g
  USE grid_dims  , ONLY : nzpmax
  USE mem_scratch, ONLY : vctr1,vctr2,vctr3,vctr4,vctr5,vctr6,vctr7, &
                          vctr8,vctr9,vctr10,vctr11,vctr12,vctr13,vctr14, &
			  vctr15,scratch
  USE mem_radiate, ONLY : ilwrtyp,iswrtyp,lonrad
  USE rconstants , ONLY : cp,cpor,p00,stefan,solar,pi180

  IMPLICIT NONE

  INTEGER :: m1,m2,m3,ia,iz,ja,jz,i,j,k,mynum

  REAL :: solfac,tdec,sdec,cdec,declin,dzsdx,dzsdy,dlon,a1,a2,dayhr,gglon  &
       ,dayhrr,hrangl,sinz,sazmut,slazim,slangl,cosi
  REAL, DIMENSION(nzpmax) :: rvr,rtr,dn0r,pird,prd,fthrl,dzmr,dztr,fthrs
  REAL, DIMENSION(nzpmax+1) :: temprd
  REAL, DIMENSION(m1,m2,m3) :: theta,pi0,pp,rv,dn0,rtp,fthrd
  REAL, DIMENSION(m2,m3) :: rtgt,f13t,f23t,glat,glon,rshort,rlong,cosz  &
       ,albedt,rlongup

 integer :: k1,k2,nlev  !-srf-rams60
  !common /radcom/ tdec,sdec,cdec,declin,rvr,rtr,dn0r,pird,prd,temprd  &
  !     ,fthrl,dzmr,dztr,fthrs

  DO j = ja,jz
     DO i = ia,iz
 !-srf-rams60
      k2=grid_g(ngrid)%lpw(i,j)
      k1=k2-1
      nlev=m1-k1+1
 !-srf-rams60

        DO k = 1,m1
!-srf-rams60
!           pird(k) = (pp(k,i,j) + pi0(k,i,j)) / cp
!	    temprd(k) = theta(k,i,j) * pird(k)
!	    rvr(k) = max(0.,rv(k,i,j))
!	    rtr(k) = max(rvr(k),rtp(k,i,j))
!	    ! Convert the next 7 variables to cgs for now.
!	    prd(k) = pird(k) ** cpor * p00 * 10.
!	    dn0r(k) = dn0(k,i,j) * 1.e-3
!	    dzmr(k) = dzm(k) / rtgt(i,j) * 1.e-2
!	    dztr(k) = dzt(k) / rtgt(i,j) * 1.e-2
!-srf-rams60
          pird(k-k1+1) = (pp(k,i,j) + pi0(k,i,j)) / cp
          temprd(k-k1+1) = theta(k,i,j) * pird(k-k1+1)
          rvr(k-k1+1) = max(0.,rv(k,i,j))
          rtr(k-k1+1) = max(rvr(k-k1+1),rtp(k,i,j))
          ! Convert the next 7 variables to cgs for now.
          prd(k-k1+1) = pird(k-k1+1) ** cpor * p00 * 10.
          dn0r(k-k1+1) = dn0(k,i,j) * 1.e-3
          dzmr(k-k1+1) = dzm(k) / rtgt(i,j) * 1.e-2
          dztr(k-k1+1) = dzt(k) / rtgt(i,j) * 1.e-2
        END DO
        temprd(1) = (rlongup(i,j) / stefan) ** 0.25
!-srf-rams60        temprd(nzp+1) = temprd(nzp)
        temprd(m1+1-k1+1) = temprd(m1-k1+1)

        ! CALL the longwave parameterizations.

        IF (ilwrtyp .eq. 2) THEN
!-srf-rams60 CALL lwradp(nzp,temprd,rvr,dn0r,dztr,pird,vctr1,fthrl,rlong(i,j))
             CALL lwradp(nlev,temprd,rvr,dn0r,dztr,pird,vctr1,fthrl,rlong(i,j))
        ELSEIF (ilwrtyp .eq. 1) THEN
!-srf-rams60 CALL lwradc(nzp+1,rvr,rtr,dn0r,temprd,prd,dztr,fthrl,rlong(i,j)  &
!-srf-rams60 	  ,vctr1,vctr2,vctr3,vctr4,vctr5,vctr6,vctr7,vctr8,vctr9,vctr10  &
!-srf-rams60 	  ,vctr11,vctr12,vctr13,vctr14,vctr15)
             CALL lwradc(nlev+1,rvr,rtr,dn0r,temprd,prd,dztr,fthrl,rlong(i,j)  &
                  ,vctr1,vctr2,vctr3,vctr4,vctr5,vctr6,vctr7,vctr8,vctr9,vctr10  &
                  ,vctr11,vctr12,vctr13,vctr14,vctr15)
        END IF

        ! The shortwave parameterizations are ONLY valid IF the cosine
        !    of the zenith angle is greater than .03 .

        IF (cosz(i,j) .gt. .03) THEN

           IF (iswrtyp .eq. 2) THEN
!-srf-rams60  CALL shradp(nzp,rvr,dn0r,dzmr,vctr1,pird,cosz(i,j)  &
!-srf-rams60 	   ,albedt(i,j),solar*1e3*solfac,fthrs,rshort(i,j))
              CALL shradp(nlev,rvr,dn0r,dzmr,vctr1,pird,cosz(i,j)  &
                   ,albedt(i,j),solar*1e3*solfac,fthrs,rshort(i,j))
           ELSEIF (iswrtyp .eq. 1) THEN


!-srf-rams60  CALL shradc(nzp+1,rvr,rtr,dn0r,dztr,prd,scratch%scr1(:)  &
!-srf-rams60       ,albedt(i,j),solar*1.e3*solfac,cosz(i,j),fthrs,rshort(i,j))
              CALL shradc(nlev+1,rvr,rtr,dn0r,dztr,prd,scratch%scr1(:)  &
                   ,albedt(i,j),solar*1.e3*solfac,cosz(i,j),fthrs,rshort(i,j))
           END IF

           ! ModIFy the DOwnward surface shortwave flux by considering
           !    the slope of the topography.

           IF (itopo .eq. 1) THEN
              dzsdx = f13t(i,j) * rtgt(i,j)
              dzsdy = f23t(i,j) * rtgt(i,j)

              ! The y- and x-directions must be true north and east for
              ! this correction. the following rotates the model y/x
              ! to the true north/east.

              ! The following rotation seems to be incorrect, so CALL this instead:
              !      SUBROUTINE uvtoueve(u,v,ue,ve,qlat,qlon,platn(ngrid),plonn(ngrid))

              dlon = (plonn(ngrid) - glon(i,j)) * pi180
              a1 = dzsdx*cos(dlon) + dzsdy * sin(dlon)
              a2 = -dzsdx*sin(dlon) + dzsdy * cos(dlon)
              dzsdx = a1
              dzsdy = a2

              dayhr = time / 3600. + float(itime1/100)  &
                   + float(mod(itime1,100)) / 60.
              gglon = glon(i,j)
              IF (lonrad .eq. 0) gglon = centlon(1)
              dayhrr = mod(dayhr+gglon/15.+24.,24.)
              hrangl = 15. * (dayhrr - 12.) * pi180
              sinz = sqrt(1. - cosz(i,j) ** 2)
              sazmut = asin(max(-1.,min(1.,cdec*sin(hrangl)/sinz)))
              IF (abs(dzsdx) .lt. 1e-20) dzsdx = 1.e-20
              IF (abs(dzsdy) .lt. 1e-20) dzsdy = 1.e-20
              slazim = 1.571 - atan2(dzsdy,dzsdx)
              slangl = atan(sqrt(dzsdx*dzsdx+dzsdy*dzsdy))
              cosi = cos(slangl) * cosz(i,j) + sin(slangl) * sinz  &
                   * cos(sazmut-slazim)
              !srf- mod rams60
	      if (cosi > 0.) then
                 rshort(i,j) = rshort(i,j) * cosi / cosz(i,j)
              else
                 rshort(i,j) =  0.
              endif
              !rshort(i,j) = rshort(i,j) * cosi / cosz(i,j)
	      !srf- mod rams60
           END IF

        ELSE
           DO k = 1,nzp
              fthrs(k) = 0.
           END DO
           rshort(i,j) = 0.
        END IF


!-srf-rams60	    DO k = 2,m1-1
!-srf-rams60	       fthrd(k,i,j) = fthrl(k) + fthrs(k)
!-srf-rams60        END DO
      ! Add fluxes, adjusting for adap if necessary.
        do k = k2,m1-1
            fthrd(k,i,j) = fthrl(k-k1+1) + fthrs(k-k1+1)
        enddo


        ! Convert the DOwnward flux at the ground to SI.

        rshort(i,j) = rshort(i,j) * 1.e-3 / (1. - albedt(i,j))
        rlong(i,j) = rlong(i,j) * 1.e-3
!-srf-rams60  fthrd(1,i,j) = fthrd(2,i,j)
      fthrd(k1,i,j) = fthrd(k2,i,j)

     END DO
  END DO
  RETURN
END SUBROUTINE radcomp

!******************************************************************************

SUBROUTINE radcomp3(m1,m2,m3,ia,iz,ja,jz,lpw  &
     ,glat,rtgt,topt,albedt,cosz,rlongup,rshort,rlong  &
     ,rv,dn0,fthrd,pi0,pp,theta,rcp)

  USE mem_grid   , ONLY: maxnzp,if_adap,zm,zt,time,ngrid
  USE mem_radiate, ONLY: iswrtyp,ilwrtyp
  USE rconstants , ONLY: cpi,p00,cpor
  USE micphys    , ONLY: level,pwmas,pwmasi,cparm,emb0,emb1,gnu,dnfac,jhcat, &
                         press,tair,emb,cx

  IMPLICIT NONE

  INTEGER :: m1,m2,m3,ia,iz,ja,jz,mcat,i,j,k
  INTEGER, DIMENSION(m2,m3) :: lpw

  REAL :: cfmasi,cparmi,glg,glgm,gammln,picpi
  REAL, DIMENSION(m2,m3) :: glat,rtgt,topt,cosz,albedt,rlongup,rshort,rlong
  REAL, DIMENSION(m1,m2,m3) :: dn0,rv,fthrd,pi0,pp,theta,rcp

  IF (level .le. 1) THEN
     mcat = 0
  ELSE
     mcat = 1
     pwmas(1) = 3.
     pwmasi(1) = 1. / pwmas(1)
     cfmasi = 1. / 524.
     cparmi = 1. / cparm
     emb0(1) = 5.24e-16
     emb1(1) = 3.35e-11
     glg = gammln(gnu(1))
     glgm = gammln(gnu(1) + pwmas(1))
     dnfac(1) = (cfmasi * exp(glg - glgm)) ** pwmasi(1)
     DO k = 2,m1-1
        jhcat(k,1) = 1
     END DO
  END IF

  DO j = ja,jz
     DO i = ia,iz

!srf-rams60        DO k = 2,m1-1
!srf-rams60           picpi = (pi0(k,i,j) + pp(k,i,j)) * cpi
!srf-rams60           press(k) = p00 * picpi ** cpor
!srf-rams60          tair(k) = theta(k,i,j) * picpi
!srf-rams60           emb(k,1) = max(emb0(1),min(emb1(1),rcp(k,i,j) * cparmi))
!srf-rams60           cx(k,1) = rcp(k,i,j) / emb(k,1)
!srf-rams60        END DO
      do k = 2,m1-1
         picpi = (pi0(k,i,j) + pp(k,i,j)) * cpi
         press(k) = p00 * picpi ** cpor
         tair(k) = theta(k,i,j) * picpi
      enddo

      if (level >= 2) then
         do k = 2,m1-1
            emb(k,1) = max(emb0(1),min(emb1(1),rcp(k,i,j) * cparmi))
            cx(k,1) = rcp(k,i,j) / emb(k,1)
         enddo
      endif

        CALL radcalc3(m1,maxnzp,mcat,iswrtyp,ilwrtyp,IF_adap,lpw(i,j)  &
             ,glat(i,j),rtgt(i,j),topt(i,j),albedt(i,j),cosz(i,j)  &
             ,rlongup(i,j),rshort(i,j),rlong(i,j)  &
             ,zm,zt,rv(1,i,j),dn0(1,i,j),fthrd(1,i,j),i,j,time,ngrid)

     END DO
  END DO
  RETURN
END SUBROUTINE radcomp3

!******************************************************************************

SUBROUTINE zen(m2,m3,ia,iz,ja,jz,jday,glat,glon,cosz,solfac)

  USE mem_grid   , ONLY: nzpmax,imonth1,idate1,iyear1,time,itime1,centlat, &
                         centlon
  USE mem_radiate, ONLY: lonrad
  USE rconstants , ONLY: pi180

  IMPLICIT NONE

  INTEGER :: m2,m3,ia,iz,ja,jz,jday,i,j,julday

  REAL :: solfac,tdec,sdec,cdec,declin,d0,d02,dayhr,radlat,cslcsd,snlsnd  &
       ,gglon,dayhrr,hrangl
  REAL, DIMENSION(nzpmax) :: rvr,rtr,dn0r,pird,prd,fthrl,dzmr,dztr,fthrs
  REAL, DIMENSION(nzpmax+1) :: temprd
  REAL, DIMENSION(m2,m3) :: glat,glon,cosz

  !common /radcom/ tdec,sdec,cdec,declin,rvr,rtr,dn0r,pird,prd,temprd  &
  !     ,fthrl,dzmr,dztr,fthrs

  jday = julday(imonth1,idate1,iyear1)
  jday = jday + nint(time/86400.)
  !      sdec - sine of declination, cdec - cosine of declination
  declin = -23.5 * cos(6.283 / 365. * (jday + 9)) * pi180
  sdec = sin(declin)
  cdec = cos(declin)

  ! Find the factor, solfac, to multiply the solar constant to correct
  ! for Earth's varying distance to the sun.

  d0 = 6.2831853 * float(jday-1) / 365.
  d02 = d0 * 2.
  solfac = 1.000110 + 0.034221 * cos (d0) + 0.001280 * sin(d0)  &
       + 0.000719 * cos(d02) + 0.000077 * sin(d02)

  ! Find the hour angle, THEN get cosine of zenith angle.

  dayhr = time / 3600. + float(itime1/100) + float(mod(itime1,100)) / 60.

  DO j = ja,jz
     DO i = ia,iz
        radlat = glat(i,j) * pi180
        IF (lonrad .eq. 0) radlat = centlat(1) * pi180
        IF (radlat .eq. declin) radlat = radlat + 1.e-5
        cslcsd = cos(radlat) * cdec
        snlsnd = sin(radlat) * sdec
        gglon = glon(i,j)
        IF (lonrad .eq. 0) gglon = centlon(1)
        dayhrr = mod(dayhr+gglon/15.+24.,24.)
        hrangl = 15. * (dayhrr - 12.) * pi180
        cosz(i,j) = snlsnd + cslcsd * cos(hrangl)
!-srf - cosz > 1 no SX6
        cosz(i,j) = min(cosz(i,j)+1.0E-10, 1.0) !LFR: Prevent 90 degrees z angle
        !cosz(i,j) = min(cosz(i,j), 1.0)
        !cosz(i,j) = max(cosz(i,j),-1.0)
!-srf	
     END DO
  END DO
  RETURN
END SUBROUTINE zen

!****************************************************************************

! SUBROUTINE radcalc3: column driver for twostream radiation code

! variables USEd within SUBROUTINE radcalc3:
! ==================================================================

! Variables in rrad3 module parameter statement

!  mb               : maximum allowed number of bands [=8]
!  mg                  : maximum allowed number of gases [=3]
!  mk               : maximum number of pseuDObands allowed for any gas [=7]
!  ncog             : number of fit coefficients (omega and asym) [=5]
!  ncb              : number of fit coefficients (extinction) [=2]
!  npartob          : number of hydrometeor categories (including dIFferent habits)
!  npartg           : number of hydrometeor categories USEd for gc coefficients [=7]
!  nrad                  : total number of radiation levels USEd (m1 - 1 + narad)
!  narad            : number of radiation levels added above model
!  nsolb            : active number of solar bands
!  nb               : active number of bands
!  ng                : active number of gases
!  jday             : julian day
!  solfac           : solar constant multiplier for variable E-S distance
!  ralcs (mb)       : rayleigh scattering integration constants
!  solar1 (mb)      : solar fluxes at top of atmosphere - corrected for ES distance
!  solar0 (mb)      : solar fluxes at top of atmosphere - uncorrected for ES distance
!  nuum (mb)        :    continuum flags
!  a0,a1,a2,a3 (mb) : Planck function fit coefficients
!  npsb (mg,mb)     : number of pseuDO bands
!  trf (mg,mb)      : reference temperature for xp and wght coefficients
!  prf (mg,mb)      : reference pressure for xp and wght coefficients
!  ulim (mg,mb)     : upper bound on pathlength for gases
!  xp (mg,mk,mb)    : coefficient USEd in computing gaseous optical depth
!  alpha (mg,mk,mb) : pressure correction factor exponent
!  beta (mg,mk,mb)  : temperature correction factor exponent
!  wght (mg,mk,mb)  : pseuDO band weight
!  exptabc (150)    : table of exponential function values
!  ocoef(ncog,mb,npartob)  : fit coefficients for hyd. single scatter.
!  bcoef(ncb,mb ,npartob)  : fit coefficients for hyd. extinction coefficient.
!  gcoef(ncog,mb,npartg)   : fit coefficients for hyd. asymmetry parameter.

! Input variables from model

!  m1               : number of vertical levels in model grid
!  ncat             : max number of hydrometeor categories [=7]
!  mcat             : actual number of hydrometeor categories [= 0, 1, or 7]
!  nhcat            : number of hydrometeor categories including ice habits [=15]
!  iswrtyp          : shortwave radiation parameterization selection flag
!  ilwrtyp          : longwave radiation parameterization selection flag
!  glat             : latitude
!  rtgt             : terrain-following coordinate metric factor
!  topt             : topography height
!  albedt          : surface albeDO
!  cosz             : solar zenith angle
!  rlongup          : upward longwave radiation at surface (W/m^2)
!  rshort           : DOwnward shortwave radiation at surface (W/m^2)
!  rlong            : DOwnward longwave radiation at surface (W/m^2)
!  jnmb (ncat)      : microphysics category flag
!  dnfac (nhcat)    : factor for computing dn from emb
!  pwmasi (nhcat)   : inverse of mass power law exponent for hydrometeors
!  zm (m1)          : model physical heights of W points (m)
!  zt (m1)          : model physical heights of T points (m)
!  press (nzpmax)   : model pressure (Pa)
!  tair (nzpmax)    : model temperature (K)
!  rv (m1)          : model vapor mixing ratio (kg/kg)
!  dn0 (m1)         : model air density (kg/m^3)
!  fthrd (m1)       : theta_il tENDency from radiation
!  jhcat (nzpmax,ncat)  : microphysics category array
!  cx (nzpmax,ncat) : hydrometeor number concentration (#/kg)
!  emb (nzpmax,ncat): hydrometeor mean mass (kg)

! Variables input from model scratch space (redefined internally on each CALL)

!  zml (nrad)       : physical heights of W points of all radiation levels (m)
!  ztl (nrad)       : physical heights of T points of all radiation levels (m)
!  dzl (nrad)       : delta-z (m) of all radiation levels
!  pl (nrad)        : pressure (Pa)
!  tl (nrad)        : temperature (K)
!  dl (nrad)        : air density of all radiation levels (kg/m^3)
!  rl (nrad)        : vapor density of all radiation levels (kg/m^3)
!  vp (nrad)        : vapor pressure (Pa)
!  o3l (nrad)       : stores the calculated ozone profile (g/m^3)
!  flxu (nrad)      : Total upwelling flux (W/m^2)
!  flxd (nrad)      : Total DOwnwelling flux (W/m^2)
!  t (nrad)         : layer transmission function
!  r (nrad)         : layer reflection function
!  tc (nrad)        : cumulative optical depth
!  sigu (nrad)      : upwelling layer source function
!  sigd (nrad)      : DOwnwelling layer source function
!  re (nrad)        : cumulative reflection function
!  vd (nrad)        : multi-scat. dIFfUSE DOwnwelling contributions
!                         from source functions
!  td (nrad)        : inverse of cumulative transmission fnct
!  vu (nrad)        : multi-scat. dIFfUSE upwelling contributions
!                         from source functions
!  tg (nrad)        : gaseous optical depth
!  tcr (nrad)       : continuum/Rayleigh optical depth
!  src (nrad)       : Planck function source for each band
!  fu (nrad*6)      : upwelling fluxes for pseuDO-bands (W/m^2)
!  fd (nrad*6)      : DOwnwelling fluxes for pseuDO-bands (W/m^2)
!  u (nrad*mg)      : path-length for gases (H_2O, CO_2, O_3)  (Pa)
!  tp (nrad*mb)     : optical depth of hydrometeors (m^-1)
!  omgp (nrad*mb)   : Single scatter albeDO of hydrometeors
!  gp (nrad*mb)     : Asymmetry factor of hydrometeors

! LoCALLy-defined variables

!  ngass (mg)       : flags indicating IF H20, CO2, O3 are active for solar wavelengths
!  ngast (mg)       : flags indicating IF H20, CO2, O3 are active for long wavelengths
!  prsnz,prsnzp     : pressure in top two model reference state levels

! Additional variables USEd ONLY within SUBROUTINE mclatchy:
! ==================================================================

!  namax            : maximum allowed number of added rad levels above model top[=10]
!                       USEd for oc and bc coefficients [=13]
! mcdat (33,9,6)    : Mclatchy sounding DATA (33 levels, 9 soundings, 6 vars)
! mclat (33,9,6)    : mcdat interpolated by season to latitude bands
! mcol (33,6)       : mclat interpolated to lat-lon of grid column

! Additional variables USEd ONLY within SUBROUTINE cloud_opt:
! ==================================================================

!  ib .......... band number
!  dn .......... characteristic diameter (m)
!  oc .......... scattering albeDO fit coefficients
!  bc .......... extinction fit coefficients
!  gc .......... asymmetery fit coefficients
!  kradcat ..... cross-reference table giving Jerry's 13 hydrometeor category
!                   numbers as a function of 15 microphysics category numbers

! Particle Numbers describe the following particles:

!     Harrington radiation code             RAMS microphysics
! ----------------------------------------------------------------
!  1:   cloud drops                 1.  cloud drops
!  2:   rain                        2.  rain
!  3:   pristine ice columns        3.  pristine ice columns
!  4:   pristine ice rosettes       4.  snow columns
!  5:   pristine ice plates         5.  aggregates
!  6:   snow columns                6.  graupel
!  7:   snow rosettes               7.  hail
!  8:   snow plates                 8.  pristine ice hexagonal plates
!  9:   aggregates columns          9.  pristine ice dENDrites
!  10:  aggregates rosettes        10.  pristine ice needles
!  11:  aggregates plates          11.  pristine ice rosettes
!  12:  graupel                    12.  snow hexagonal plates
!  13:  hail                       13.  snow dENDrites
!                                  14.  snow needles
!                                  15.  snow rosettes

! for the asymmetery parameter, since we ONLY have spherical
! particles, there are ONLY 7 particle types...

!  1:   cloud drops
!  2:   rain
!  3:   pristine ice
!  4:   snow
!  5:   aggregates
!  6:   graupel
!  7:   hail

!******************************************************************************

SUBROUTINE radcalc3(m1,maxnzp,mcat,iswrtyp,ilwrtyp,IF_adap,lpw  &
     ,glat,rtgt,topt,albedt,cosz,rlongup,rshort,rlong  &
     ,zm,zt,rv,dn0,fthrd,i,j,time,ngrid)

  USE rconstants, ONLY: cp
  USE rrad3     , ONLY: mg,namax,mb,nrad,narad,nb,ocoef,bcoef,gcoef,ncog,ncb, &
                        npartob,npartg,ng,nsolb,npsb,xp,alpha,beta,wght,prf, &
		        trf,ralcs,solar1,ulim,nuum,a0,a1,a2,a3,exptabc
  USE micphys   , ONLY: press,tair

  IMPLICIT NONE

  INTEGER :: m1,maxnzp,mcat,ilwrtyp,iswrtyp,IF_adap,lpw,ngrid
  INTEGER :: i,j,k,kk,koff,nradmax  !,ik,ib,ig

  ! ORIGINAL
  !INTEGER, SAVE :: nCALL = 0
  ! ModIFied by Alvaro L.FazENDa:
  ! SAVE in line 734 automatiCALLy SAVE all variables
  INTEGER :: ncall = 0

  INTEGER ngass(mg),ngast(mg)
  REAL prsnz,prsnzp

  REAL glat,rtgt,topt,cosz,albedt,rlongup,rshort,rlong,slr,time  !,rmix
  REAL eps   !,dzl9,rvk1,rvk0

  REAL zm(m1),zt(m1),dn0(m1),rv(m1),fthrd(m1)

  REAL, allocatable, DIMENSION(:)     :: zml,ztl,dzl,pl,tl,dl,rl,o3l  &
       ,vp,flxu,flxd,tg,tcr,src,t,r,tc  &
       ,sigu,sigd,re,vd,td,vu  &
       ,u,fu,fd,tp,omgp,gp

  DATA eps/1.e-15/

  !     one can choose the gases of importance here,
  !       ngas = 1    gas active
  !            = 0    gas not active
  !
  !       ngas(1) = H2O
  !       ngas(2) = CO2
  !       ngas(3) =  O3

  DATA ngass/1, 1, 1/,ngast/1, 1, 1/
  SAVE

  IF (nCALL == 0) THEN
     nCALL = 1
     nradmax = maxnzp + namax
     ALLOCATE(zml(nradmax),ztl (nradmax),dzl (nradmax),pl (nradmax)  &
          ,tl (nradmax),dl  (nradmax),rl  (nradmax),o3l(nradmax)  &
          ,vp (nradmax),flxu(nradmax),flxd(nradmax),tg (nradmax)  &
          ,tcr(nradmax),src (nradmax),t   (nradmax),r  (nradmax)  &
          ,tc (nradmax),sigu(nradmax),sigd(nradmax),re (nradmax)  &
          ,vd (nradmax),td  (nradmax),vu  (nradmax)               )
     ALLOCATE(u(nradmax*mg),fu(nradmax*6),fd(nradmax*6))
     ALLOCATE(tp(nradmax*mb),omgp(nradmax*mb),gp(nradmax*mb))

     !**ALF: Putting zero on all variables ALLOCATEd:
     CALL azero5(nradmax, zml, ztl, dzl, pl, tl)
     CALL azero5(nradmax, dl, rl, o3l, vp, flxu)
     CALL azero5(nradmax, flxd, tg, tcr, src, t)
     CALL azero5(nradmax, r, tc, sigu, sigd, re)
     CALL azero3(nradmax, vd, td, vu)

     CALL azero(nradmax*mg, u)

     CALL azero2(nradmax*6, fu, fd)

     CALL azero3(nradmax*mb, tp, omgp, gp)

  END IF


  koff = lpw - 2
  nrad = m1 - 1 + narad - koff

  CALL mclatchy(3,m1,IF_adap,koff  &
       ,prsnz,prsnzp,glat,rtgt,topt,rlongup  &
       ,zm,zt,press,tair,dn0,rv,zml,ztl,pl,tl,dl,rl,o3l,dzl)

  ! zero out scratch arrays

  CALL azero(nrad*mg,u)
  CALL azero(nrad*6,fu)
  CALL azero(nrad*6,fd)
  CALL azero(nrad*mb,tp)
  CALL azero(nrad*mb,omgp)
  CALL azero(nrad*mb,gp)

  ! Compute hydrometeor optical properties

  CALL cloud_opt(mb,nb,nrad,m1,koff,mcat,dzl  &
       ,dn0,tp,omgp,gp &
       ,ocoef,bcoef,gcoef,ncog,ncb,npartob,npartg,i,j,time)

  ! Get the path lengths for the various gases...

  CALL path_lengths(nrad,u,rl,dzl,dl,o3l,vp,pl,eps)

  DO k = 1,nrad
     IF (rl(k) .lt.   0. .or.  &
          dl(k) .lt.   0. .or.  &
          pl(k) .lt.   0. .or.  &
          o3l(k) .lt.   0. .or.  &
          tl(k) .lt. 183.) THEN

        PRINT*, 'Temperature too low or negative value of'
        PRINT*, 'density, vapor, pressure, or ozone'
        PRINT*, 'when CALLing Harrington radiation'
        PRINT*, 'at k,i,j = ',k,i,j,'   ngrid=',ngrid
        PRINT*, 'STOPping model'
        PRINT*, 'rad: rl(k), dl(k), pl(k), o3l(k), tl(k)'
        DO kk=1,nrad
           WRITE (*,'(5g15.6)') rl(kk), dl(kk), pl(kk), o3l(kk), tl(kk)
        END DO
        STOP 'STOP: radiation CALL'
     END IF
  END DO

  ! CALL shortwave and longwave schemes...

  IF (iswrtyp .eq. 3 .and. cosz .gt. 0.03) THEN
     CALL azero2(nrad,flxu,flxd)

     CALL swrad(nrad,ng,nb,nsolb,npsb,   &         !  counters
          u,pl,tl,dzl,vp,                  &      !  model variables
          xp,alpha,beta,wght,prf,trf,ralcs,  &   !  band specIFics
          solar1,ngass,                      &    !        "
          albedt,slr,cosz,               & !  boundaries
          tg,tcr,tp,omgp,gp,   &!  loCALLy defined
          t,r,tc,                              & !  for fluxes
          sigu,sigd,re,vd,td,vu,fu,fd,         & !       "
          flxu,flxd,ulim,i,time)                     !  sw fluxes

     rshort = flxd(1)
!srf-rams60     DO k = lpw,m1-1
!srf-rams60        fthrd(k+koff) = fthrd(k+koff)  &
!srf-rams60             + (flxd(k) - flxd(k-1) + flxu(k-1) - flxu(k)) / (dl(k) * dzl(k) * cp)
!srf-rams60     END DO
     do k = lpw,m1-1
       kk = k - koff
       fthrd(k) = fthrd(k)  &
          + (flxd(kk) - flxd(kk-1) + flxu(kk-1) - flxu(kk)) &
             / (dl(kk) * dzl(kk) * cp)
     enddo
  ELSE
     rshort = 0.
  END IF

  IF (ilwrtyp .eq. 3) THEN
     CALL azero2(nrad,flxu,flxd)

     CALL lwrad(nrad,ng,nb,nsolb,npsb,nuum,   &    !  counters
          u,pl,tl,dzl,vp,                       & !  model variables
          xp,alpha,beta,wght,prf,trf,ralcs,     &!  band specIFics
          a0,a1,a2,a3,                          &!        "
          exptabc,ngast,                        &!  boundaries
          tg,tcr,tp,omgp,gp,src,   &  !  loCALLy defined
          t,r,tc,                               &!  for fluxes
          sigu,sigd,re,vd,td,vu,fu,fd,          &!       "
          flxu,flxd,ulim)                          !  fluxes, output

     rlong = flxd(1)
!srf-rams60 	 DO k = lpw,m1-1
!srf-rams60 	    fthrd(k+koff) = fthrd(k+koff)  &
!srf-rams60 		 + (flxd(k) - flxd(k-1) + flxu(k-1) - flxu(k)) / (dl(k) * dzl(k) * cp)
!srf-rams60     END DO
   do k = lpw,m1-1
      kk = k - koff
      fthrd(k) = fthrd(k)  &
         + (flxd(kk) - flxd(kk-1) + flxu(kk-1) - flxu(kk)) &
            / (dl(kk) * dzl(kk) * cp)
   enddo

  END IF

  RETURN
END SUBROUTINE radcalc3

!******************************************************************************

SUBROUTINE cloud_opt(mb,nb,nrad,m1,koff,mcat,dzl  &
     ,dn0,tp,omgp,gp,oc,bc,gc,ncog,ncb,npartob,npartg,i,j,time)

  ! computing properties of spherical liquid water and irregular ice
  ! using fits to adt theory
  !
  ! ib .......... band number
  ! mb .......... maximum number of bands
  ! nb .......... total number of bands
  ! m1 .......... number of vertical levels
  ! dzl .......... delta z in each level (m)
  ! dn .......... characteristic diameter (m)
  ! emb ......... mean hydrometeor mass (kg)
  ! cx .......... hydrometeor concentration (#/kg)
  ! tp .......... optical depth
  ! omgp ........ scattering albeDO
  ! gp .......... asymmetry parameter
  ! oc .......... scattering albeDO fit coefficients
  ! bc .......... extinction fit coefficients
  ! gc .......... asymmetry fit coefficients
  ! ncog ........ number of fit coefficients (omega and asym)
  ! ncb ......... number of fit coefficients (extinction)
  ! kradcat ..... cross-reference table giving Jerry's 13 hydrometeor category
  !                 numbers as a function of 15 microphysics category numbers
  ! dn0 ......... model air density (kg/m^3)
  ! dnfac ....... factor for computing dn from emb
  ! pwmasi ...... inverse of power USEd in mass power law
  ! npartob ..... number of hydrometeor categories (including dIFferent habits)
  !                 USEd for oc and bc coefficients
  ! npartg ...... number of hydrometeor categories USEd for gc coefficients
  !

  USE micphys, ONLY: jnmb,cx,jhcat,dnfac,emb,pwmasi

  IMPLICIT NONE

  INTEGER mb,nb,ib,nrad,m1,ncog,ncb,krc,npartob,npartg,koff  !,iz
  INTEGER icat,mcat,k,i,j,ihcat

  INTEGER kradcat(15)
  REAL dzl(nrad),tp(nrad,mb),omgp(nrad,mb),gp(nrad,mb),dn0(m1)
  REAL oc(ncog,mb,npartob),bc(ncb,mb,npartob),gc(ncog,mb,npartg)
  REAL ext,om,gg,dn,time

  REAL dnmin(7),dnmax(7)
  DATA dnmin /   1.,   10.,   1.,  125.,   10.,   10.,   10./
  DATA dnmax /1000.,10000., 125.,10000.,10000.,10000.,10000./

  DATA kradcat/1,2,3,6,10,12,13,5,5,3,4,8,8,6,7/

  DO icat = 1,mcat
     IF (jnmb(icat) .gt. 0) THEN

        DO k = 2,m1-1-koff

           IF (cx(k+koff,icat) .gt. 1.e-9) THEN
              ihcat = jhcat(k+koff,icat)
              krc = kradcat(ihcat)
              dn = dnfac(ihcat) * emb(k+koff,icat) ** pwmasi(ihcat) * 1.e6
              dn = max(dnmin(icat),min(dnmax(icat),dn))

              DO ib = 1,nb

                 ext = cx(k+koff,icat) * dn0(k+koff) * dzl(k)  &
                      * bc(1,ib,krc) * dn ** bc(2,ib,krc)
                 om = oc(1,ib,krc)  &
                      + oc(2,ib,krc) * exp(oc(3,ib,krc) * dn)  &
                      + oc(4,ib,krc) * exp(oc(5,ib,krc) * dn)
                 gg = gc(1,ib,icat)  &
                      + gc(2,ib,icat) * exp(gc(3,ib,icat) * dn)  &
                      + gc(4,ib,icat) * exp(gc(5,ib,icat) * dn)


                 tp(k,ib) = tp(k,ib) + ext

                 omgp(k,ib) = omgp(k,ib) + om * ext
                 gp(k,ib) = gp(k,ib) + gg * om * ext

              END DO

           END IF
        END DO

     END IF
  END DO

  ! Combine the optical properties....

  DO ib = 1,nb
     DO k = 2,m1-1-koff
        IF (tp(k,ib) .gt. 0.0) THEN
           gp(k,ib) = gp(k,ib) / omgp(k,ib)
           omgp(k,ib) = omgp(k,ib) / tp(k,ib)
        ELSE
           omgp(k,ib) = 0.0
           gp(k,ib) = 0.0
        END IF

        ! Check for validity of opt values before CALLing radiation
        !      IF (tp(k,ib) .lt. 0) THEN
        !         PRINT*, 'tp(k,ib) less than zero for k,ib = ',k,ib
        !         PRINT*, 'tp(k,ib) = ',tp(k,ib)
        !         STOP 'opt1'
        !      END IF
        !      IF (omgp(k,ib) .lt. 0. .or. omgp(k,ib) .gt. 1.) THEN
        !         PRINT*, 'omgp(k,ib) out of range [0,1] for k,ib = ',k,ib
        !         PRINT*, 'omgp(k,ib) = ',omgp(k,ib)
        !         STOP 'opt2'
        !      END IF
        !      IF (gp(k,ib) .lt. 0. .or. gp(k,ib) .gt. 1.) THEN
        !         PRINT*, 'gp(k,ib) out of range [0,1] for k,ib = ',k,ib
        !         PRINT*, 'gp(k,ib) = ',gp(k,ib)
        !         STOP 'opt3'
        !      END IF

     END DO
  END DO

  RETURN
END SUBROUTINE cloud_opt

!------------------------------------------------------------------

SUBROUTINE path_lengths(nrad,u,rl,dzl,dl,o3l,vp,pl,eps)

  ! Get the path lengths for the various gases...

  IMPLICIT NONE
  INTEGER :: nrad
  REAL, DIMENSION(nrad) :: rl,dzl,dl,o3l,vp,pl
  REAL, DIMENSION(nrad,3) :: u
  REAL :: rvk0,rvk1,dzl9,rmix,eps
  INTEGER :: k

  u(1,1) = .5 * (rl(2) + rl(1)) * 9.81 * dzl(1)
  u(1,2) = .5 * (dl(2) + dl(1)) * .45575e-3 * 9.81 * dzl(1)
  u(1,3) = o3l(1) * 9.81 * dzl(1)

  rvk0 = rl(1)
  DO k = 2,nrad
     rvk1 = (rl(k) + 1.e-6)
     dzl9 = 9.81 * dzl(k)
     rmix = rvk1 / dl(k)
     vp(k) = pl(k) * rmix / (.622 + rmix)
     u(k,1) = (rvk1 - rvk0) / (log(rvk1 / rvk0) + eps) * dzl9
     u(k,2) = (dl(k) - dl(k-1)) / (log(dl(k) / dl(k-1)) + eps)  &
          * dzl9 * 0.45575e-3
     u(k,3) = 0.5 * dzl9 * (o3l(k) + o3l(k-1))
     rvk0 = rvk1
  END DO

  RETURN
!:DOC%END
END SUBROUTINE path_lengths




